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ABSTRACT 

In this paper we describe an algorithm and deduce the related mathematical formulae 
that allows the computation of observed fluxes in stellar and planetary systems with 
arbitrary number of bodies being part of a transit or occultation event. The presented 
method does not have any limits or constraints for the geometry and can be applied 
for almost all of the available limb darkening models as well. As a demonstration, 
we apply this scheme to gather information for the orbital inclinations from multiple 
transiting planetary systems in cases when mutual transits occur. We also show here 
that these mutual events constrain inclinations unambiguously, yielding a complete 
picture for the whole system. 

Key words: Celestial mechanics - Stars: Binaries: Eclipsing - Stars: Planetary Sys- 
tems - Methods: Analytical - Techniques: Photometric 



1 INTRODUCTION 

Wit h the adv ent of the space missions Corot and Kepler 
fsee lBarge et al. 2008; Borucki et al. 2009) and as the result 
of numerous successful ground-based surveys, nearly two- 
hundred transiting extrasolar planets are known up to date. 
Furthermore, the number of candidates awaiting for confirm- 
ing t he planetary propert ies exceeds the magnitude of thou- 
sand (|Borucki et al.ll201ll ). These systems gives us a unique 
perspective for various studies because most of the planetary 
and orbital parameters can be obtained without any ambi- 
guity. Transiting planet ary companions a re also known in 
multiple stellar systems jPovle et al.l l201l'). Moreover, both 
systems of transiting planets and eclipsing binaries provide 
substantial information about the stars from which the ab- 
solute physical properties can be easily obtained, completely 
independently from other methods and therefore these stud- 
ies are essential to confirm stellar evolution models. 

In the case of multiple transiting planetary systems (see 
e.g. lHolmarill2010l ). triple or hier archical stellar sys tems or 
circumbinary planetary systems (|Dovle et al.ll201ll ). plane- 
tary systems around one of a binary component s or systems 
with planetary companions (exomoons , see e.g. ISzabo et al.l 
l2006l : 1Simon et al.ll2009l : lKippin!dl2009l ). there is a chance to 
observe mutual eclipsing or transiting events when (at least) 
three of the bodies are aligned along the line of sight. More- 
over, if planet and/or stellar formation prefers co- planar or- 
bits, the chance is even higher. Conversely, observing mutual 
transit events yields additional information ab out the orbital 
characteristics of the whole system. Recentlv. lSato fc Asadal 
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f200<f) and 'Ragozzine fc HolmanI (|2010l ) analyzed these ef- 
fects and their qualitative influence on extrasolar moons and 
multiple planetary systems. In this paper we discuss how the 
photometric measurements are affected due to such mutual 
transiting or eclipsing events by giving an algorithm that 
models light curves of such phenomena. Our method de- 
scribed here is capable to compute light curve models for 
arbitrary number of eclipsing or transiting bodies and for 
all of the well-known limb darkening models without any 
restrictions for the projected diameters of the active compo- 
nents. 

The structure of this paper is as follows. Section [5] de- 
scribes the algorithm and the formulae needed to evaluate 
the fluxes or light curve points for events with multiple tran- 
siting, eclipsing or occulting companions. In Section [3] we 
briefly discuss the qualitative properties of systems where 
mutual transits may occur and demonstrate how informa- 
tion gathered from such mutual events can be exploited in 
order to constrain orbital alignments in such a transiting 
planetary system. And finally, Section[3]summarizes the key 
points and results of this paper. 



2 THE LIGHT CURVE MODEL 

In this section we briefly describe the methods used to com- 
pute the light curve models for multiple transiting objects. 
Recently, Kipping (20_f I) published an algorithm that is ca- 
pable to estimate the observed flux when two bodies transit 
their host star simultaneously. However, that method works 
only when one of these bodies is very small (i.e. assuming 
a homogeneous flux density beyond this very small disk). 
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Here we demonstrate an alternative algorithm that is signif- 
icantly more con cise and can be treated as an extension of 
the approach by iKippind (|201lh in several ways. First, the 
presented method is capable to incorporate more than two 
transiting bodies. Space-borne missions l ike Kepler are ex- 
pected t o detect both extrasolar moons (jSzabo et al.ll2006l : 
iKippingj I2OO9I ) by different methods (e.g. detecting timing 
variations or via photometry) and systems with three or 
more transitin g planetary compani ons that are also known 
(Kepler-11, see Lissauer et aLllioill ). Second, the model can 
be extended for various limb darkening models, that can be 
quantified by a series expansion on the apparent stellar sur- 
face (and might lack circular symmetry). Such models can 
also be exploited to q uantify cases with asymmetric light 
curves like K0f-13(b) (jSzabo et al.ll201ll l. Third, in terms 
of computation time, the method presented here can also be 
an alternative for the well-known models available for the 
single -planet cases (jMandel fc Agoi2002l : [GimeneJ2006l : [pi] 
I2OO8I ). Also, this method do es not require several dozens of 
distinct geometric cases (see iMandel fc A"goill2002l : [ Kipping 
l20lj ). And finally, more sophisticated cases like non-uniform 
thermal radiations can also be considered, even in mutually 
transiting systems. This might be relevant in the analysis of 
near-infrared light curves of close-in eclipsing companions. 

The computation of the presented light model is based 
on three subsequent steps. First, a net of disjoint arcs is ob- 
tained from the mutual intersections of the apparent stellar 
and/or planetary disc edges. Second, we generate a vector 
field whose exterior derivative (i.e. the planar component 
of the curl operator) is the surface brightness. The surface 
brightness must be in accordance with our assumptions for 
the limb darkening models. Third, we apply Green's or the 
Kelvin-Stokes theorem (known from differential geometry or 
vector calculus) to integrate the vector field on an appropri- 
ately oriented subset the arcs. In the following, we discuss 
these three steps as well as their applications for various 
surface brightness functions. 



2.1 Net of arcs 

In principle, the projected stellar, planetary or lunar discs 
are characterized by the center coordinates xo,yo and the 
radius r. An arc on one of these circles is quantified by the 
additional parameters ip'^'^^ and Aip, where V.(«' is the posi- 
tion angle between the reference axis (a;-|-) and the beginning 
of the arc while is the length of the arc in radians. All of 
the arcs in this model are oriented in counter-clockwise (i.e. 
prograde or positive) direction. In the following, the circles 
and arcs are indexed by k and I, respectively. Obviously, for 
the fcth circle, 



and 

(0) 



2-K 



(0) 



k,t- 



(1) 



(2) 



Completely disjoint circles or circles of which edge does not 
intersect other ones have only one arc that is used to rep- 
resent the circle itself. Namely, {£} = {1} and /S.ipk,i ~ 27r 
while the value of ip^^\ can be arbitrary. 

The net of arcs is built iteratively. If the new circle is 
disjoint, only one arc is placed with Aipk, 1 = 27r, otherwise 




Figure 1. A complex configuration of 5 circles. The arcs are 
denoted by k:{Cy (} where k and k' are the indices of the circles 
and the set Cy ^ is a list of circle indices in which the given arc 
goes through. The thick arcs mark the boundary of the region 
in circle #1 that is disjoint from the other circles. It is easy to 
see that these boundary arcs are labelled by either 1:{} or 
(where 2 ^ fc). 



the two position angles for the two intersection points are 
computed using known trigonometric relations and the ap- 
propriate arcs are split into two or three smaller ones. After 
obtaining this set of arcs, the topology is also generated. 
Namely, by checking for each arc what are the circles that 
contains this arc inside. Let us denote this subset of circles 
regarding to the ^th arc of the fcth circle by Ck,e- Note that 
this set might be empty or might even contain all of the 
circle indices with the exception of fc. See e.g. Fig. [T] for a 
particular example of 5 intersecting circles. 

2.2 Surface brightness and vector fields 

The surface brightness of the host star can be modelled 
by various limb darkening lawfl Recalling Green's theorem 



we can write 



(D Af)dA 



f • dr. 



(3) 



as 



DAf = 



(4) 



Here S C R and dS is the boundary of S. These are two 
and one-dimensional manifolds on on which the standard 
measures are dA and dr, respectively. This equation can also 
be viewed as a somewhat special case to the Stokes theorem 
known for the curl operator in a three-dimensional space. 
The term D A f denotes the exterior derivative of f , that 
can be written in terms of vector components as 

dfy df^ 
dx dy 

For our problem discussed in this paper we can apply the 
above equation © as follows. First, we have to find a func- 
tion f = {fx,fy) of which exterior derivative is the given 
stellar surface brightness density. We have to note that due 
to the Young theorem, this function is ambiguous, since we 
can add an arbitrary scalar gradient to f , of which addition 

bee e.g. ' http:// www . astro, keele. ac . uk/j kt /codes /j kt Id . html | 
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does not change its exterior derivative. Therefore, it is rec- 
ommended to find such an f which have "nice properties" 
making the computation of the integral on the right-hand 
side of equation (|3]) convenient. For instance, a homogeneous 
surface can be modelled with the function 



fi 



1 / -V 

2\+x 



(5) 



The exterior derivative of this function is unity and there is 
no preferred direction or position angle in the vector field 
described by fi . 



2.3 Integration on the arcs 

Let us consider a set arcs of which union is the boundary: 
a C dS. The arc a corresponding to the circle centered at 
{xa,ya) with a radius of ra and parameterized via the posi- 
tion angle (p implies the measure 



dr= I -"^^"^|d^. 

-r cos (p 



(6) 



Thus, the total fiux F coming from the area S is then com- 
puted as 



F : 



E 

aeas 



[fy (x, y) cos ip- f^{x,y) sin (p] r^dt/p, 



(7) 



where for more compact notations, we define x = Xa + 
Ta cos ip and y = ya + ra sin pi. Note that the integration 
limits (pi"'' and (p^a^ are not necessarily p>'^^ and p'^a^ + Apa, 
because the direction of the ^ f -dr integral must be positive 
in all cases. If we use (pi^^ — p'^^ and p>^a^ — ip^a^ + Atpa, we 
have to multiply the integrand by ±1, depending whether 
the right-hand directed arc a points inside or outside the 
area S (see e.g. Fig. [1] or Fig. [2] for examples and further 
explanation). 



2.4 Some surface brightness functions 

Now, we compute the integrals behind the sum of equa- 
tion ([7]) for various surface brightness functions. Let us con- 
sider a star whose projected disk center is located at (0, 0) 
and has a radius of unity. The domain of the functions of 
out interest is this unit circle, i.e. 4- 1. 



2.4-1 Homogeneous surface 

As we have seen earlier (see e.g. equation [5]), the vector field 
f = {—\y, has a curl of unity. For a given arc a = 0, 

the integrals behind the sum of equation ((?]) can be written 



Fq 



i(2;o + r cos (p)r cos p + ^(yo + r sin (p)r sin p ■ 



V2 



1 , . \ 1 2 

-r(xocosip + yosmp) + -r 



r{p2 



P2) + -rxo{smp2 



sin^i) -I- 



+ ^ryo{cospi 



cos pi2 ) 



(8) 



Note that the value of Fq does depend on the actual choice 
for f , i.e. it would be different if we add a gradient to the 
vector field f . However, J]] Fa in equation ((Tjl would not be 

a 

altered after such an addition of a gradient field. By sum- 
ming the results yielded by equation ((S]) for the arcs {a}, 
we can easil y rep r oduce the results of in Sections 2, 3 and 
Fig. 5 



1 easil y repi r oduce 
of lKippind l|201lh . 



2.4-2 Polynomial intensities 

Various limb darkening models contain terms which can be 
quantified as polynomial functions of the (a;, y) centroid co- 
ordinates (for instance, the quadratic limb darkening law). 
In additional, any analytical limb darkening profiles can be 
well approximated by polynomial functions, therefore it is 
worth to compute terms in equation (O for such cases. 

Without any restrictions, let us consider the term x^y"^ . 
Due to the linearity of the integral and summation, if the 
surface intensity can be described by polynomials, comput- 
ing the integral in equation (O for the above terms are suf- 
ficient. First, let us define 



Mr, 



(xq + r cos p>)^{yo + r sin p>Y 



(9) 



and introduce c = cos 93 and a — sin(/p, just for simplicity. 
Thus, in the expansion of equation we should compute 
expressions like J MpqC or J MpqS. Here we give a set of 
recurrence relations with which these indefinite integrals can 
be evaluated. It is easy to show that 

/ Mpq = xoj Mp-i,q +rj Mp_i,,c, or (10) 

J Mpq = yo J Mp,q-l + r J Mp^q^lS. (11) 

For the terms J MpqC and J MpqS we can write 

{I + p + q) j' MpqC = +MpqS + rpJ Mp-i,q+ (12) 
+XOP J Mp-i,qC + yoqJ Mp,q_ic, 
-MpqC + rqj Mp,,_i -f- (13) 
+XOP J Mp-i,qS + yoq J Mp^q-is. 



{1+p + q) j MpqS 



In order to bootstrap these set of recurrence relations, we 
only have to use the following: 



Moo = 1, 

/Moo = /l = id, 

/ MqoC = +M00S, 

J Moos = -Mooc. 



(14) 
(15) 
(16) 
(17) 



However, for some cases we might compute these inte- 
grals a bit more easier. For instance, the surface density 
x^ + y^ can be integrated as the exterior derivative of 
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Figure 2. A "movie" of a mutual transit, caused by two relatively large companion. The host star is the largest fixed circle while the two 
companions are the smaller circles moving from up to down. The arc boundaries of the area(s) on which the surface intensity is integrated 
are denoted by thick lines. Thick solid lines mark where the integrals in equation (O are directed in prograde (counter-clockwise) direction 
while thick dashed lines mark where the integrals are directed in retrograde (clockwise) direction. Thin lines mark the other arcs, which 
are irrelevant regarding to the integration. Note that even in the topologically complex cases, the number of relevant boundary arcs is 
three or less (with the exception of the sixth frame where the number of relevant arcs is six) . 



f = i^—^x^y — ij/^j+ix^ + ^xy^Y Hence, the primitive in- 
tegral in equation ((7|) for this f will be 

F[>f] = {24(xg + yo')r + 12r=')^ - (18) 

~4yo{6xl + 2yl + 9r^) cos<^ + 
+Axo{2xi + Qyt-, + <dr) sm if — 
— 4r^[j/o cos(3v5) + a;o sin(3<^)] — 
— 24a;oyo'~ cos(2(^) — r^sin(4((9)} . 



2.4.3 Linear hmb darkening 

The well known formula of linear limb darkening law gives 
us a surface flux density that can be written in the form of 

I{x,y) ^1- 41 ~ 11], (19) 

where /i = y^l — x^ — y^ and c denotes the linear limb dark- 
ening coefficient. Performing integrals on this function will 
yield a linear combination of integrating a constant (see 
Section 12.4. ip and integrating the function y^l — x'^ ~ y^. 
Thus, here we derive a function of which exterior derivative 
is y'l — x^ — y^ as well as we compute the indefinite inte- 
gral that is going to appear in equation ((7)|. First of all, we 
have to note that this function is defined only in the domain 
bounded by the circle x^ + y'^ — 1. 

Let us search the function f of which exterior deriva- 
tive is ^1 — x^ — y'^ = \/l — r^ in the form of f = 
[~yf(r^)/'2,+xfir'^)/2]. For simplicity, we introduce r = 
x'^ + y'^. It is easy to show that the expansion of equa- 
tion (|4)| for this function yields an ordinary differential equa- 
tion (ODE) for /(■) that is 



/(0+C,f'(C) = Vl^, (20) 
where ^ = r^ . One of the solutions of this ODE is 

This solution behaves analytically at ^ = 0, i.e. at the center 
of the stellar disk. Therefore, the function f can be written 
as 



1-(1- 



3(a;2 -f y2) 



+x 



(22) 



By substituting this x — {fx, fy) into equation (0, one finds 
that the primitive integral 



1 /■ 1 - (1 - - 2rpcosip 
3 J R'^ + 2rp cos (p' 



-(r 4- pcos (p')rd(f' 



(23) 

should be computed, if we parameterize the arc (on which 
the f is integrated) similarly as earlier. The constants are 
the following: = Xq + yg, E? = r^ + p^, and (p' is defined 
to be p cos (fi' equals to xq cos ip + yo sin ip. The primitive in- 
tegral in equation 1)23^ can be computed analytically and 
this computation yields a function which contains elemen- 
tary functions as well as elliptic integrals. As a first step, let 
us introduce the constants 



92 

S2 

d2 

Q 

SE 

riE 

F 
E 
P 



r^ + p^ + 2rp cos (p' , 

(r + pf, 

(r- pf, 
1 

y/rp' 



2cos{p'/2)^ 



rp 



da 



2V rp ' 

1 - d2 

d2 ' 

F(se; fefi), 
E(se; fefi), 
n(sE; riE, /ce) 



(24) 
(25) 
(26) 

(27) 
(28) 

(29) 

(30) 

(31) 
(32) 
(33) 



Here F(-;-), E(-;-) and n(-; ■, •) denotes the incomplete ellip- 
tic integrals of the first, second and third kind, respectively. 
Then, F[ip'] in equation (|23|l is computed as 



F[p'] 



■ arctan 



2 2 
p — r 



+ 



+ — + -rp^T~q^smp' + 
5 9 

+ i(l-4r^ + 2/)QF + 

+ irp(7r^ + 5rp + p^ -4)QF + 

+ ^rp{8-Ur^ -2p^)QE + 



+ 



Ir + p 
6 r — p 



QP. 



(34) 



One may note s ome similarities between t hese terms and 
the equations in iMandel fc Agoll (2002) or^M (|2008l '). Al- 
though the formulae above include incomplete elliptic inte- 
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Figure 3. Simulated light curves of transits caused by two, relatively large planets and occurring nearly simultaneously but vfithout any 
overlap (left panel). A light curve of a mutual transit with almost the same geometry is displayed on the right panel. The signal- to- noise 
ratio on these plots are nearly the same as one expects from Kepler photometry. The thick solid lines show the expected flux i/ the 
two simultaneous transits would be treated independently, i.e. the apparent overlapping of the planets would be neglected. See text for 
further details, including the geometric parameters of these simulated transits. 



grals, the actual evaluation of these does not require longer 
computation time than the complete ones. Both types of 
elliptic integrals are computed vi a the Carlson symmetric 
forms (|Carlson fc Gustafson|[l993l ). for which computation 
very fast and rob ust algorithms are available in the litera- 
ture (|Press et al.| [i992: Carlson 1994). 

It should also be noted that the evaluation of equa- 
tion (|34p might be done with caution in some cases where 
the values of the variables or constants defined in equa- 
tions (I24|l - (|33|) introduce singularities in some of the terms. 
These values correspond to cases where the arc endpoints are 
at the edge of the bounding circle at x^+y^ = 1 and/or when 
the arc intersect the origin (i.e. \ix = y = Q). However, these 
singularities yield more simple formulae in general. For in- 
stance, the case oi p = Q (i.e. when the bounding circle and 
the arc is concentric), equation (I34|l becomes simply 

F[^'] = i[l-(l-r^)^/^]^'. (35) 

Especially, when the arc is a part of the bounding circle itself 
(which is a practically frequent case: see e.g. the thick solid 
lines in the frames of Fig. [2]), i.e. when r = 1 and p = 0, 
then -F [</?'] is merely <y5'/3- these cases should be 

treated carefully during a practical implementation. 



2.4.^ Quadratic limb darkening 

The quadratic limb darkening stellar profile is characterized 
by the surface flux density 7 = l — ci{l — fj,) — C2(l — fi) . Since 
fi — y^l — — y^, by expanding this equation, we obtain a 
constant term, with a value of 1 — Ci — 2c2, a polynomial term 
x'^ + y^ with a coefficient +C2 and a term that is proportional 
to fi and has a coefficient ci -|- 2c2. Hence, the formulae in 
the previous three subsections (|2.4.1lE".4.2l and l2.4.3p can be 
applied accordingly to evaluate the final apparent fiuxes in 
the case of a quadratic limb darkening law. 



3 ORBITAL INCLINATIONS 

Mutual transits occur when at least two bodies (that can 
be, for instance, two planets or a planet and its moon) tran- 
sits the host star simultaneously and their projections also 
overlap. Due to this overlap, the observed flux coming from 
the host star is larger than if we would consider naively the 
flux decreases from each body independently. In Fig. O we 
display a series of images that clearly show this effect. In the 
previous section we deduced the algorithms and mathemat- 
ical formulae that are needed for the computation of the 
total observed flux for arbitrary geometry and for various 
limb darkening models. 

As a demonstration, in Fig. [3] we display two sim- 
ulated light curves with nearly the same orbital geome- 
try. The planet-to-size ratio for the two companions are 
Ri/R-, = 0.13 and R2/Ri, = 0.10 while the orbital pa- 
rameters are the following: a\/Ri, = 4.3000, 61 = 0.35, 
nx = 2.0d-\ aa/i?,, = 9.5952, 62 = 0.22, na = O.Gd'S 
ASl = 18° and both of the planets have a circular orbit. Here 
rifc denotes the orbital angular frequency: it is = 2-K/Pk, 
where Pk is the orbital period, bk is the impact parame- 
ter of the transit, ak/R* is the normalized semimajor axis 
(in the units of stellar radii) and ASl = Q.i — 0,2 is differ- 
ence in the orbital ascending nodes (note that the reference 
plane here is the plane of the sky) . The mid-transit time of 
the inner planet is Ei = 0.02 d while the outer planet has 
E2 = —0.06 d on the left panel, and E2 — 0.00 d on the right 
panel. This difference between the mid-transit times yields a 
mutual transit in the latter case (see the flux excess in Fig. [3] 
aX t ~ 0.02 . . . 0.08 d) while there is no overlap between the 
apparent planetary disks in the former case. 

It can easily be seen that the time evolution of the flux 
excess yielded by the mutual transi10 has similar qualitative 
properties as the normal transits have. Namely, it has a mid- 



^ Here we treat this "flux excess" relative to the flux level that 
would be if we neglect the effect of the overlapping and simply 
calculate the yield of the two components independently. 
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time, a peak and a duration. The larger the flux excess peak, 
the larger the overlapping area is. At a first glance, the only 
quantity for which an observation of a mutual transit yields 
additional constraints is the difference in the Af2, the dif- 
ference between the orbital ascending nodes. Qualitatively, 
the longer the duration of this flux excess, the smaller the 
absolute value of Af2 i^. However, the depth and the ex- 
act time of the mutual event defines the impact parameters 
more precisely. This is rather relevant when one or both of 
the impact parameters are relatively small: the uncertainty 
of iP does not strongly depend o n the actual value of h (see 
e.g. Pal 2008:: ICarter et aklliooi ). thus the uncertainty in 6 
will be rather large for small values of h due to the relation 
A6 = (26)~^A(fe^). Indeed, for instance, the analysis of the 
light curves shown in Fig. [3] yields the following. If no mu- 
tual transit occurs (left panel), the best-fit values for fc^'s 
wiU be &i = 0.323 ±0.012 and 62 = 0.215 ±0.016 while if we 
can observe the mutual transit, we obtain 61 = 0.351 ±0.003 
and 62 = 0.223 ± 0.007 while for the node difference we got 
AJ7 = 17.4±0.5°. For this demonstration of light curve anal- 
ysis, we employed an improved Markov Chain Monte-Carlo 
algorithm as implemented in the If it utility (Pal 2009). 

Of course, if the difference in the nodes, Af2 = Q.i — 0.2 
is known, we can compute the mutual inclination im of the 
orbits as well using the well-known relation 



= cos i\ cos 12 ± sin i\ sin 12 cos Af2. 



(36) 



It should also be mentioned that the analysis of mutual tran- 
sits resolve the ambiguity between the values of ±Af2. And 
of course, the precise analysis of mutual transits should in- 
volve the gravitati onal interactions between the companions 
(see e.g. lpSllioiol ). especially when data are available on a 
timescale on which the perturbations are not negligible (con- 
trary to the demonstration presented here). 



4 DISCUSSION 

In this paper we investigated the possibilities for comput- 
ing apparent stellar fluxes in multiple or hierarchical stellar 
and/or planetary systems during simultaneous transits or 
occultations. The presented algorithm is capable to derive 
these fluxes for arbitrary number of bodies that are actively 
parts of the transiting or eclipsing event. This method can 
then be applied for various analyses of complex astrophysi- 
cal systems, including multiple transiting planetary systems, 
hierarchical stellar systems with planetary companions and 
extrasolar moons as well. 

Currently, the algorithm is implemented in 
ANSI C, in the form of a plug-in module for 
the program If it and available from the address 
http : //szof i . elte.hu/~apal/utils/astro/mttr/. This 
module features functions named mttrXyC), where X 
denotes the number of transiting bodies and y can be 
"u", "1" or "q" for the uniform flux density, linear limb 
darkening and quadratic limb darkening. Evidently, these 
functions have 3X ± y parameters where j/ is 0, 1 or 2 
for "u", "1" or "q", respectively. The current version of 

^ Imagine two completely retrograde orbits: in this case, the rel- 
ative speed of the transiting planets is the highest, thus the du- 
ration of the overlapping will be the smallest. 



this module does not compute the parametric derivatives 
of the functions analytically but emulates them using 
numerical approximations for the If it utility. Since both 
the parametric derivatives of the arcs (with respect to the 
circle center coordinates and radii) and the parametric 
derivatives of equation ((T)) can be computed analytically, 
the composition of these two would give us the required 
derivatives. 

As a demonstration, we applied this method to obtain 
mutual inclinations of orbits in multiple transiting planetary 
systems. The analysis presented here clearly shows that ob- 
serving a mutual transits yields not only an accurate value 
for the ascending node difference but also results a more 
precise value for the impact parameters, and therefore the 
orbital inclinations as well. 
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